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The tools of statistical mechanics and nonlinear dy- 
namics have been used successfully in the past to ana- 
lyze models of social phenomena ranging from language 
choice [ill to political party affiliation |2] to war [2| and 
peace In this work, we focus on social systems com- 
prised of two mutually exclusive groups in competition 
for members [5HlO) . We compile and analyze a new data 
set quantifying the declining rates of religious affiliation 
in a variety of regions worldwide and present a theory to 
explain this trend. 

People claiming no religious affiliation constitute the 
fastest growing religious minority in many countries 
throughout the world [TT]. Americans without religious 
affiliation comprise the only religious group growing in 
all 50 states; in 2008 those claiming no religion rose to 
15 percent nationwide, with a maximum in Vermont at 
34 percent 12 . In the Netherlands nearly half the popu- 
lation is religiously unaffiliated [13]. Here we use a min- 
imal model of competition for members between social 
groups to explain historical census data on the growth of 
religious non-affiliation in 85 regions around the world. 
According to the model, a single parameter quantifying 
the perceived utility of adhering to a religion determines 
whether the unaffiliated group will grow in a society. The 
model predicts that for societies in which the perceived 
utility of not adhering is greater than the utility of ad- 
hering, religion will be driven toward extinction. 

MODEL 

We begin by idealizing a society as partitioned into 
two mutually exclusive social groups, X and Y, the un- 
affiliated and those who adhere to a religion. We as- 
sume the attractiveness of a group increases with the 



number of members, which is consistent with research 
on social conformity [T3HI3- We further assume that at- 
tractiveness also increases with the perceived utility of 
the group, a quantity encompassing many factors includ- 
ing the social, economic, political and security benefits 
derived from membership as well as spiritual or moral 
consonance with a group. Then a simple model of the 
dynamics of conversion is given by[l] 

— = yPyx{x,u^) - xPxy{x,Ux) (1) 

where Pyx{x, u^) is the probability, per unit of time, that 
an individual converts from Y to X, x is the fraction of 
the population adhering to X at time t, < m^; < 1 
is a measure of AT's perceived utility, and y and Uy 
are complementary fractions to x and u^. We require 
) — Pyx{^ — x,l — Ux) to obtain symmetry un- 
der exchange of x and y and Pyx(x,0) = Pyx{0,Ux) = 
to capture the idea that no one will switch to a group 
with no utility or adherents. The assumptions regarding 
the attractiveness of a social group also imply that Pyx is 
smooth and monotonically increasing in both arguments. 
Under these assumptions, for generic Pyx{x,Ux) Eq. ([T]) 
has at most three fixed points, the stability of which de- 
pends on the details of Pyx (see Supplementary Material 
Section S2). 

A functional form for the transition probabilities con- 
sistent with the minimal assumptions of the model is 
Pyx{x,Ux) = cx°'Ux, where c and a are constants that 
scale time and determine the relative importance of x 
and Ux in attracting converts, respectively. (Supplemen- 
tary Figure S2 illustrates the structure of the fixed points 
for this case.) If a > 1 there are three fixed points, one 
each at X = and x = 1, which are stable, and one at 
< X < 1, which is unstable. For a < 1 the stability of 
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these fixed points is reversed. For the boundary at a = 1, 
there are only two fixed points, one of which is stable and 
the other unstable (see Supplementary Material Section 
S3). 
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FIG. 1. Percentage religiously unaffiliated versus time in four 
regions: (a) the autonomous Aland islands region of Finland, 
(b) Schwyz Canton in Switzerland, (c) Vienna Province in 
Austria, (d) the Netherlands. Red dots indicate data points 
from census surveys, black lines indicate model fits. Relative 
utilities for the religiously unaffiliated populations as deter- 
mined by model fits were Ux — 0.63, 0.70, 0.58, 0.56. 

In Figure [T] we fit the model to historical census data 
from regions of Finland, Switzerland, Austria, and the 
Netherlands, four of 85 worldwide locations for which we 
compiled and analyzed data. The initial fraction unaf- 
filiated Xq and the perceived utility were varied to 
optimize the fit to each data set, while c and a were 
taken to be global. A broad minimum in the error near 
a = 1 indicated that as a reasonable choice (see Supple- 
mentary Material Section S4). Figure [T](d) shows that, 
if the model is accurate, nearly 70% of the Netherlands 
will be non-affiliated by midcentury. 

The behavior of the model can be understood ana- 
lytically for a = 1, in which case we have dx/dt = 
cx{l — x){2ux — 1): logistic growth. An analysis of the 
fixed points of this equation tells us that religion will 
disappear if its perceived utility is less than that of non- 
affiliation, regardless of how large a fraction initially ad- 
heres to a religion. However, if a is less than but close 
to one, a small social group can indefinitely coexist with 
a large social group. Even if a > 1 it is possible that so- 
ciety will reach such a state if model assumptions break 
down when the population is nearly all one group. 

Figure [2] shows the totality of the data collected and 
a comparison to the prediction of Eq. ([I]) with a — 1, 
demonstrating the general agreement with our model. 
Time has been rescaled in each data set and the ori- 
gin shifted so that they lie on top of one another. See 
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FIG. 2. All data on changes in religious affiliation with time 
(85 data sets). Time has been rescaled so data sets lie on top 
of one another and the solution curve with Ux — 0.65. Red 
dots correspond to regions within countries, while blue dots 
correspond to entire countries. Black line indicates model 
prediction for Ux = 0.65. 

Supplementary Material Section S3 for more details. 

Our assumption that the perceived utility of a social 
group remains constant may be approximately true for 
long stretches of time, but there may also be abrupt 
changes in perceived utility, a possibility that is not in- 
cluded in the model. We speculate that for most of hu- 
man history, the perceived utility of religion was high 
and of non- affiliation low. Religiously non-affiliated peo- 
ple persisted but in small numbers. With the birth of 
modern secular societies, the perceived utility of adher- 
ence to religion versus non-affiliation has changed signifi- 
cantly in numerous countries^Il , such as those with cen- 
sus data shown in Fig. [ij and the United States, where 
non-affiliation is growing rapidly ^lE]. 

One might ask whether our model explains data better 
than a simple empirical curve. Logistic growth would be 
a reasonable null hypothesis for the observed data, but 
here we have provided a theoretical framework for ex- 
pecting a more general growth law ([T]), and have shown 
that data suggest logistic growth as a particular case of 
the general law. Our framework includes a rational math- 
ematical foundation for the observed growth law. 

GENERALIZATIONS 

We have thus far assumed that society is highly inter- 
connected in the sense that individual benefits stem from 
membership in the group that has an overall majority. 
For that reason, the model as written is best applied on 
a small spatial scale where interaction is more nearly all- 
to-all. We can generalize this model to include the effects 
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of social networks: rather than an individual deriving 
benefits from membership in the global majority group, 
he or she will instead benefit from belonging to the local 
majority among his or her social contacts[71 [12]. In or- 
der to define "local," however, we must introduce either 
a spatial dimension to the problem, or a network defining 
social interaction. On a network, Eq. ([!]) becomes: 

d{R, 



dt 
where 



= (1 - {Ri))Pyx{Xi,Ux) - {Ri)Pyx{i - Xi,l- Ux) 

(2) 



N IN 

x,= ^A,,i?, /^A,j (3) 

i=i / i=i 

defines the local mean religious affiliation, A is a binary 
adjacency matrix defining the social network and i? is a 
binary religious affiliation vector (1 indicates membership 
in group X). An ensemble average has been assumed in 
order to write a derivative for the expected religious affil- 
iation {Ri{t)) in ([2]), since this system is stochastic rather 
than deterministic. In the all-to-all coupling limit, A = 1 
and X 2 — X ^ SO ^ reduces to Q. 

A further generalization to a continuous system with 
arbitrary coupling can be constructed with the introduc- 
tion of a spatial dimension. The spatial coordinate ^ will 
be allowed to vary from —1 to 1 with a normalized cou- 
pling kernel G(^, C) determining the strength of connec- 
tion between spatial coordinates ^ and The religious 
affiliation variable R now varies spatially and temporally 
with < R{^,t) < 1, so individuals may have varying 
degrees of affiliation. Then the dynamics of R satisfy 

dR 

— = (1 - R)Pyxix, Ux) - RPyxil -X,l- Ux) (4) 

in analogy with the discrete system. Here x represents 
the local mean religious affiliation. 



x{^,t) 



G{^,e)R{e,m' 



(5) 



Note that simulation of Eq. ^ with continuous real- 
valued R and large N is equivalent to integration of 
Eq. Q with appropriate initial conditions and appro- 
priately chosen G(^,^'). This is because ([2| goes from a 
stochastic system for binary R, to a deterministic system 
for real i? G [0, 1]. 

In the case of all-to-all coupling, G(^,^') = 1/2, and 
x{£_,t) = i f_iR{Cit)d^' = Rit), independent of space, 
where R is the spatially averag ed value of R. Then Q 
becomes 

dR 

— = (1 - R)PyxiR, Ux) - RPyxil ~R,l~Ux) . (6) 

If at some time t R{^,t) = i?o(i) is independent of 
space, then R{t) = Ro{t) and Eq. ^ becomes 

OR 

= il-Ro)Pyx{Ro,Ux)-RoPyA^-RoA-Ux) , (7) 



which follows dynamics identical to the original two- 
group discrete system ([T]). 

We can impose perturbations to both the coupling 
kernel (i.e., the social network structure) and the spa- 
tial distribution of R values to examine the stability of 
this system and the robustness of our results for the all- 
to-all case. One very destabilizing perturbation consists 
of perturbing the system towards two separate clusters 
with different R values. These clusters might represent a 
polarized society that consists of two social cliques in 
which members of each clique are more strongly con- 
nected to members of their clique than to members of 
the other clique. Mathematically, this can be written as 
G(^,e') = i + i<5(2i/(0 - l)(2i/(e') - 1), where <5 is 
a small parameter ((5 <C 1) that determines the ampli- 
tude of the perturbation and H{£^) represents the Heav- 
iside step function. This kernel implies that individuals 
with the same sign of ^ are more strongly coupled to one- 
another than they are to individuals with opposite-signed 

The above perturbation alone is not sufficient to 
change the dynamics of the system — a uniform state 
P{^yto) = Ro will still evolve according to the dynamics 
of the original system ([T]) (See Supplementary Material 
Section S5). 

We add a further perturbation to the spatial distri- 
bution by imposing i?(^,to) = Ro + ^ sgn(^), where e is 
a small parameter. This should conspire with the per- 
turbed coupling kernel to maximally destabilize the uni- 
form state. 

Surprisingly, an analysis of the resulting dynamics re- 
veals that this perturbed system must ultimately tend to 
the same steady state as the original unperturbed sys- 
tem. Furthermore, the spatial perturbation must even- 
tually decay exponentially, although an initial growth is 
possible. (See Supplementary Material Section S5 for 
more details on the perturbative analysis.) 

The implication of this analysis is that systems that are 
nearly all-to-all should behave very similarly to an all-to- 
all system. In the next Section we describe a numerical 
experiment that tests this prediction. 



NUMERICAL EXPERIMENT 

We design our experiment with the goal of controlling 
the perturbation from an all-to-all network through a sin- 
gle parameter. We construct a social network consisting 
of two all-to-all clusters initially disconnected from one- 
another, and then add links between any two nodes in 
opposite clusters with probability p. Thus p — I cor- 
responds to an all-to-all network that should simulate 
([!]), while p = leaves the network with two discon- 
nected components. Small perturbations from all-to-all 
correspond to p near 1, and p can be related to the cou- 
pling kernel perturbation parameter S described above as 
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FIG. 3. Results of simulation of the discrete stochastic model 
([2| on a network with two initial clusters weakly coupled to 
one another. The ratio p of out-group coupling strength to 
in-group coupling strength is (a) p — 0.01; (b) p — 0.40; (c) 
p — 0.80 (A^ = 500). Steady states are nearly identical to the 
predictions of the all-to-all model ([TJ. 
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FIG. 4. Variation in the behavior of system Q with increas- 
ing perturbation off of all-to-all (A*' — 500, a;o = 0.1, m = 0.6). 
Equivalent values of the perturbation parameter 5 in order of 
decreasing p are S = 0, S = 0.14, 5 — 0.60, and 5 — 0.98. 



p = {l — 6)/{l + 5) (assuming all links in the network have 
equal weight). The size of each cluster is determined by 
the initial condition xq as Nx — xqN, Ny = (1 — xo)N. 

Figure [3] compares the results of simulation of system 
([2]) with varying perturbations off of all-to-all. The theo- 
retical (all-to-all) separatrix between basins of attraction 
is a vertical line at — 1/2. Even when p = 0.01, when 
in-group connections are 100 times more numerous than 
out-group connections, the steady states of the system 
and basins of attraction remain essentially unchanged. 

In the case of the continuous deterministic system Q , 
the equivalent figure to [3] is extremely boring: numeri- 
cally, the steady states of the perturbed system are in- 
distinguishable from those of the unperturbed all-to-all 
system, regardless of the value of p (see Supplementary 
Section S6 and Figure S5). 

The only notable difference between the dynamics of 



the continuous networked system and the dynamics of 
the original all-to-all system ([T]) is a time delay d appar- 
ent before the onset of significant shift between groups 
(see Figure |4|. We were able to find an approximate ex- 
pression for that time delay as d esc — \np/{2u — 1) (see 
Supplementary Material Section S7, figures S6 and S7). 

What we have shown by the generalization of the 
model to include network structure is surprising: even 
if conformity to a local majority influences group mem- 
bership, the existence of some out-group connections is 
enough to drive one group to dominance and the other to 
extinction. In the language of references the pop- 

ulation will reach the same consensus, despite the exis- 
tence of individual cliques, as it would without cliques, 
with only the addition of a time delay. 

In a modern secular society there are many opportuni- 
ties for out-group connections to form due to the preva- 
lence of socially integrated institutions — schools, work- 
places, recreational clubs, etc. Our analysis shows that 
just a few out-group connections are sufficient to explain 
the good fit of Eq. ([T]) to data, even though Eq. ([!]) im- 
plicitly assumes all-to-all coupling. 



CONCLUSIONS 

We have developed a general framework for model- 
ing competition between social groups and analyzed the 
behavior of the model under modest assumptions. We 
found that a particular case of the solution fits census 
data on competition between religious and irreligious seg- 
ments of modern secular societies in 85 regions around 
the world. The model indicates that in these societies 
the perceived utility of religious non-affiliation is greater 
than that of adhering to a religion, and therefore pre- 
dicts continued growth of non-affiliation, tending toward 
the disappearance of religion. According to our calcu- 
lations, the steady-state predictions should remain valid 
under small perturbations to the all-to-all network struc- 
ture that the model assumes, and, in fact, the all-to- 
all analysis remains applicable to networks very differ- 
ent from all-to-all. Even an idealized highly polarized 
society with a two-clique network structure follows the 
dynamics of our all-to-all model closely, albeit with the 
introduction of a time delay. This perturbation analysis 
suggests why the simple all-to-all model fits data from 
societies that undoubtedly have more complex network 
structures. 

For decades, authors have commented on the surpris- 
ingly rapid decline of organized religion in many regions 
of the world. The work we have presented does not ex- 
clude previous models, but provides a new framework for 
the understanding of different models of human behavior 
in majority/minority social systems in which groups com- 
pete for members. We believe that, with the application 
of techniques from the mathematics of dynamical sys- 



terns and perturbation theory, we have gained a deeper 
understanding of how various assumptions about human 
behavior will play out in the real world. 

This work was funded by Northwestern University and 
The James S. McDonnell Foundation. The authors thank 
P. Zuckerman for useful correspondence. 
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1. GENERALITY OF MODEL 

The model presented in this paper is applied to the 
widespread phenomenon of religious shift, but may be 
more generally applicable to a variety of competitive so- 
cial systems. The model allows for either competitive ex- 
clusion (a > 1) or stable coexistence (a < 1) in systems 
composed of two social groups, and makes sense in the 
context of social networks. A similar model (reference 
10) was applied to the phenomenon of language death. 
Some other competitive social systems in which identical 
or very similar models may apply include, for example, 
smoker vs. non-smoker, vegetarian vs. meat-eater, obese 
vs. non-obese, and Mac user vs. PC user. 



2. WHY THREE FIXED POINTS 

In the main text of our paper, we state that there 
can be at most three fixed points for "generic" functions 
Pyx{x;Ux) that satisfy our assumptions of symmetry, 
monotonicity, C°° continuity, and limiting properties. In 
this section we will clarify the meaning of "generic" . 

Pyx is a non-negative function of x parametrized by Ux 
(which will henceforth be abbreviated as simply u). The 
fixed points x* can be written as solutions to the equation 
= (1 — x)Pyx(x; u) — xPyx{y — x;l — u) for a given value 



of u. When u = the limiting properties Pyx{Q]u) = 
and Pyx{x;Q) — 0, along with monotonicity, imply that 
only X* = Q and x* = \ can be fixed points, with x* = 
the only stable fixed point. When u = 1, similarly, only 
X* — Q and a;* = 1 can be fixed points, with only x* — \ 
stable. 

If there is a single intermediate fixed point x* ^ Q,x* ^ 
1 for all values of u G (0, 1), then it must limit to x* — > 
when u — >■ and a;* — >■ 1 when u — > 1 (assuming it's 
stable — the opposite will be true if it is unstable). In or- 
der for other fixed points to appear, the continuous curve 
connecting {x; u) — (0; 0) to {x; u) = (1; 1) would have to 
have zero slope at some value of u (see Supporting Figure 
SI). Thus the condition for a single intermediate fixed 
point is that dx/du > for all u (stable), or dx/du < 
for all u (unstable). 




SUPPORTING FIG. SI. Typical fixed points for Eq. (1). 
Here Pyx{x;u) — cx°'Ux, with (a) a = 3 and (b) a = |. Red 
open lines indicate unstable branches, black solid fines indi- 
cate stable branches of fixed points. Panel (a) demonstrates 
that the intermediate unstable branch of fixed points x^(ux) 
serves as a separatrix, with all other initial conditions leading 
to a; = or a; = 1. Panel (b) demonstrates how the stable 
fixed point x*s{ux) typically varies with u^. If the interme- 
diate fixed point is unstable, it must limit to a;Jj — >■ 1 when 
Ux Q and x* — >■ when — > 1. 



For a separable function Pyx{x;u) = X{x)U{u), 
X{x) > 0, U{u) > 0, the implications of this 
condition are as follows. The fixed point equation 

(1 — x)Pyx{x;u) = xPyx{l — x;l — u) becomes (1 — 
x)X{x)U{u) — xX{l — x)U{l — u), and, assuming x 
and X ^ 1, 



Thus 



U{u 


) 


X X{l-x) 


U{1- 


u) 1 


- X X{x) ■ 


\ U{u) 1 


' du 


x A(l-x)1' 


_C/(1 - u) 


dx 


l-x X{x) 



Since 



Uj-u 



u(i-u)\ = [U'{u)U{l^u)+U{u)U'{l-u)]/U\l- 
u) > V u due to assumptions (monotonicity implies that 
U and U' must both be positive for all nonzero argu- 
ments), the condition dx/du < becomes 



X{1 - x) 



X{x) 



< . 
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This can be simplified to 

X'(x) X'(l -x) 1 1 

— — H ^ > - H . 

X{x) X{l~x) X l-x 

The direction of the inequaUty would be reversed for a 
stable intermediate fixed point. Note that a sufficient 
condition would be X'/X > 1/x. This, or the same 
condition with the inequality reversed, is clearly satis- 
fied for any power law form Pyx{x;u) ~ x"", a ^ 1. 
It is also satisfied for any function with a monotonia 
first derivative X'{x) (Sketch of proof: Let X'{x) — 
Xq + f{x), where f{x) = Jq X"{£^)d^ is a monotoni- 
cally increasing function. Then xX'{x) = xXq + xf{x) 
and Xix) = J^X'iOd^ = xX^ + J^fiOd^. Thus 
xX'{x) - X{x) = xf{x) - Jo fiOd^- That last quan- 
tity is necessarily greater than zero for any monotoni- 
cally increasing f{x), and therefore xX'{x) > X{x), or 
X'/X > 1/x.) 

An analogous result holds for inseparable functions 
Pyx{x]u). Using the same approach, the condition is: 

-Pyx{x; u) + —^P'yxil -x]l~u)> 

X " I — X " 

■^Pyx{x; u) + Jl—7/^^Py^i'^ -x]l-u) , 

where prime notation represents a derivative with respect 
to the argument, not the parameter. Thus a sufficient 
condition is Pyx{x;u)/Pyx{x;u) > 1/x for all u. The full 
condition is satisfied by any function for which curvature 
doesn't change sign. 

3. STABILITY OF FIXED POINTS 

Examine the stability of the fixed point at a; = (and 
note that the same argument will work for the stability of 
the fixed point at a: = 1). Set x = rj, a small perturbation 
from X = 0. Then 

77 = (1 - ■n)Pyx{m u) - vPyxil -rr,i- u) 

« Pyx{0; u) + T] [P^,(0; u) - PyxiO; u) -Pyx{l;l- u)] 
= -,; [P,,(l;l-7.)-P;,(0;u)] 

to Oijf'). So the fixed point = is stable to small 
perturbations if Pyx{^] 1 — u) > P'yx{^]u). For a power 
law Pyx ^ x"" , this will be true only when a > 1. The 
fixed point = will be unstable for a < 0, and its 
stability will depend on the sign of u — ^ when a = 1. 

The stability of the intermediate fixed point is fully de- 
termined once the stabilities of the two endpoints a;* = 
and X* = 1 are known. Because it is a one- dimensional 
flow, the intermediate fixed point must be stable when 
the endpoints are unstable, and vice- versa when the end- 
points are stable (see Figure S2). 



a>1 ©-eO > ? >^ ? • 



a = 1 O > ^ >^ ^ > • 



a<1 O > >^ ^ •-eO 

SUPPORTING FIG. S2. The flow in x for various values 
of the constant a. Fifled circles indicate stable fixed points, 
while open circles indicate unstable fixed points. The leftmost 
fixed points correspond to s = 0, while the rightmost fixed 
points correspond to 2; = 1. 

4. DATA SETS AND MODEL FITS 

Data used in validating this model originated in census 
surveys from a range of countries worldwide. A total of 85 
data sets had 5 or more independent data points. These 
came from various regions of 9 different countries: Aus- 
tralia, Austria, Canada, the Czech Republic, Finland, 
Ireland, the Netherlands, New Zealand, and Switzerland. 

Fitting was done by minimizing root-mean-square er- 
rors. Using a functional form Pyx ~ cx°'u, the parameters 
c and a were taken to be universal while the parameter u 
was allowed to vary with each data set. This was accom- 
plished by simultaneously optimizing c, a, and ui...un 
such that the RMS error summed over all N data sets 
was minimized. 

Supporting Figure S3 shows how that summed error 
varied with the parameter a. We chose a 1 for the fits 
discussed in this paper both for simplicity and because of 
the broad minimum visible around a — 1. The parameter 
c, which simply sets a time scale, was approximately 0.2. 

5. PERTURBATION OF NETWORK 
STRUCTURE 

In this section we examine in greater depth the im- 
plications of Eq. (4), a continuous deterministic system 
with arbitrary coupling. 

All-to-all coupling 

If G(^, £,') = 1/2 then there is uniform all-to-all coupling, 
and we see that x{£,,t) — \ J^i R{Ct)d^' — Rit), inde- 
pendent of space, where R is the spatially averaged value 
olR. 



3 



4.5 



0} 

I— 

as 

IT 

03 
I 

C 

CO 
0) 

E 

J, 3-5 

O 

O 



filter used to obtain 
smoothed curve 




0.8 



SUPPORTING FIG. S3. Summed root-mean-square error 
over all data sets versus parameter a in Py^ = cx°'Ux- The er- 
ror was calculated by finding the combination of parameters c, 
xoi and u^i (where i varies over all data sets) that minimized 
the root mean square error between the model predictions and 
the data. Blue curve indicates exact error calculations, red 
indicates smoothed error after convolution with a Gaussian 
(inset). Note that there appears to be a broad minimum near 
a = 1. 



Then (4) becomes 
dR 



dt 



= {l-R)Py,{R-u)~RPy^{l-R-l~u) (SI) 



If at some time t* R{^,t*) = Ro{t*) is independent of 
space, then R{t*) = Ro{t*) and Eq. (SI) becomes 



dRo 
dt 



= (l- Ro)PyM;u) - RoPyA^- Ro-A-u) , (S2) 



which follows dynamics identical to the original two- 
group discrete system. 



Perturbation off of uniform R with all-to-all coupling 



We impose a destabilizing perturbation such that the 
portion of the population with ^ < has lower R and 
the portion with ^ > has higher R, i.e., 

i?(^,to) = i?o + esgn(e) , (S3) 

where e is a small parameter. Then x{S.) — 
\ \\ R{£.')dC = Ro, and from (4) we get 



OR 



-(l-i?o-esgn(0)P^,(i?o;u) 



-{Ro + esgn{0)Py.{l-Ro;l-u) (S4) 



We can also look at the dynamics of the mean religious 
affiliation R, 



dt dt\2 



dt 



Plugging (S4) into (S5) and simplifying gives 
dR dRa 



(S5) 



dt 



dt 



il-Ro)Py^ (Ro; u) ~ RaPy^il -Ro;l-u) , 

(S6) 

so the mean religious affiliation R{t) continues to follow 
the dynamics of the original system despite the pertur- 
bation. 

Rearranging (S4), we see 



dR dRo 



- e sgniOiPy^iR; u) + Py^^l -R;l-u)) , 

(S7) 



dt dt 

and direct differentiation of (S3) yields 
dR dRo de 



dt 



dt dt 



sgn(0 



(S8) 



Equating these expressions yields a differential equation 
for e{t): 



de 
dt 



= -e{Py,{R- u) + Py,{l ~R-A~ u)) 



(89) 



Note that e remains independent of the spatial coordi- 
nate, and that e ^ as t — > oo, for any initial condi- 
tion (the time constant may vary with the parameter u 
and the state R). So the initial perturbation must damp 
out, and the system must evolve to a single affiliation as 
t oo, just as the original system (1) did. 



Non-uniform coupling 



We consider the case of non-uniform spatial coupling as 
the continuum limit of a discrete network where the links 
are nearly but not quite all-to-all. In that case, a very 
destabilizing perturbation would be one in which the 
network is segregated into two clusters, each one more 
strongly coupled internally than externally. As described 
in the main text, one kernel representing such coupling 
is 

G(C, e) = l + l^HiO - l){2H{a ~ 1) , (SIO) 

where (5 is a small parameter {S ^ I) that determines 
the amplitude of the perturbation and -ff (^) represents 
the Heaviside step function. 

If the initial state of the population is uniform such 
that R{£.,ta) = Rq then x{£„to) = jliG{^,C)Rod£.' = 
Ro and R will satisfy Eq. (4), giving 



dRo 
dt 



{l-Ro)PyARo;u)-RoPy^il-Ro;l-u) . (Sll) 
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Thus R will remain uniform in space and will follow the 
same dynamics as the original system, despite a non- 
uniform coupling kernel of arbitrary amplitude. 



Perturbation off of uniform R with non-uniform 
coupling 



As before, we impose a destabilizing perturbation such 
that the portion of the population with ^ < has lower R 
and the portion with ^ > has higher R, i.e., R{S,,to) = 
i?o + e sgn(^), where e is again a small parameter. This 
should conspire with the perturbed coupling kernel to 
maximally destabilize the uniform state. 
Now (5) gives 



x{^,t)=iRo-e) I ^G{i,f)di'+ 



and from (4) we get 



(Ro + e) f G{U')d£.' 
Jo 

=Ro + e5 sgn(^) , 



OR 



= {l-Ro-e sgn(0)P„x(i?o + eS sgn(0;u) 



dt " ^ "^"vs/;^ yxK 

-{Ro + e sgn{0)PyM -Ro-eS sgn(0; 1-u) . 

(S13) 

Directly differentiating i?(^,to) = i?o + e sgn(^), then 
rearranging terms on the right-hand-side of (S13) gives 

+ sgn(0^ = (1 - iio)P.x(iio + e6 sgn(0; u) 

- -Ro-Pycc(l - Ro-eS sgn(^); 1 - u) - e sgn(^) [Pyx{Rf) 
+ eS sgn(^); u) + Pyx{i - Ro-eS sgn(^); 1-u)] . 

(S14) 

Now calculate ^ directly. Note that Ro = R = 
'^J^,RiC,m',so 



(S12) 



dRo _d fl 
~W ~dt V 2 



/>(c',.K')4/_; 



1 dR{£,',t)^, 



dt 



= 2 y J(l - -Ro - e sgn{0)Pyx{Ro + e6 sgn(e'); -{Ro + e sgn{0)Pyx{l -Ro-eS sgn(f ); 1 - u)\ 
= ^(1 - Ro) [Py^{Ro - eS; u) + Py^{Ro + eS; u)] - ^i?o [Pyx{l - Ro + eS;l - u) + Py^{l -Ro-eS;l- u)] 
+ ^e[Pyx(i?o - e^; u) - Py^{Ro + eS; u)] + ^e[Py^{l - Ro + eS;l - u) - Py^{l -Ro-eS;l- u)] . 



r 



Taylor expanding in both e and S eliminates all first order 
terms in e, leaving 



dRo 
dt 



= {l-Ro)Py^{Ro;u)-RoPy^{l-Ro;l-u)+0{e^S) , 

(S15) 

so the assumption that the mean religious affiliation fol- 
lows the dynamics of the original unperturbed system is 
well justified. 

Similarly Taylor expanding Eq. (S14) to first order in 
e and S allows canceling of the dRo/dt terms on either 
side (using (S15)), leaving an equation in e: 

1^ = - e^Py^{Ro; u) + Py^{l - i?o; 1 - u) 
-S[{l-Ro)P'yARa;u) 
+ RoP'y,{l-Ro;l-u)]Y (S16) 

The sign of the quantity in braces in (S16) determines 



the stability of the uniform spatial state. It's clear that 
for sufficiently small S, the uniform state will always be 
stable. However, in systems with an unstable interme- 
diate fixed point (or no intermediate fixed point), the 
uniform state will remain stable even when the quantity 
in braces is initially positive! This is because Eq. (S15) 
will still hold for small e, making Ro approach a steady 
state value Rq = Q ov Rq = 1 from the original system. 
Since e < min(i?o, 1 — Ro) is required to maintain vari- 
ables in the allowed domain, e may initially grow, but it 
will have to eventually shrink as Rq ^ Rq- 

The above further shows that e does not develop any 
additional spatial structure, so an initial state with R = 
Ro + e sgn(^) will maintain such a spatial structure as Rq 
and e evolve in time. 

This calculation demonstrates that an understanding 
of the simple all-to-all discrete system gives insight into 
the more complex problem of religious shift on a social 



5 



network. In numerical experiment, the results of the per- 
turbation calculation described here remain valid even for 
very sparse networks quite different from all-to-all. 

6. NUMERICS 

In the main text, we describe a numerical experiment 
that we performed on systems (2) and (4). Figure 3 of 
the main text shows the results of that experiment with 
a simulated size N = 500, and in Figure S4, we show 
that the all-to-all system (1) becomes a good match to 
the discrete stochastic system (2) as the number of nodes 
increases (thus explaining why Figure 3 is well predicted 
by understanding the all-to-all system). 

1| • — ^ ^ • ~^ 1 
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SUPPORTING FIG. S5. Results of simulation of the contin- 
uous deterministic system (4) on a network with two initial 
clusters weakly coupled to one another. The ratio p of out- 
group coupling strength to in-group coupling strength is (a) 
p = 0.01; (b) p = 0.40; (c) p = 0.80 (iV = 500). When 
Ux = 1/2, all points are fixed points, so the initial condition 
determines the final state. Steady states are indistinguishable 
from those of the all-to-all model (1) despite the non-uniform 
coupling and inhomogeneous initial conditions. 
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SUPPORTING FIG. S4. Comparison of simulation of discrete 
stochastic system (2) to model predictions for various system 
sizes (all-to-all coupling). Here xq = 0.1 and the total size of 
the network is 50 (red), 100 (blue), or 500 (green). The solid 
black line represents the solution to Eq. (1), the large A'' limit 
of this model. 

We also performed a our numerical experiment on (4), 
the continuous deterministic generalization of (1). Re- 
sults are presented in Figure S5, where the steady states 
are indistinguishable from the all-to-all model. 



time tc when a solution R{tc) = \{\ -\- Rq) and tcg when 
Raii-to-au{tco) = ^i^ + Ro)- We observe this quantity to 
increase monotonically with the perturbation off of all-to- 
all S — see Figure S6 for typical behavior at various values 
of S. In the limit that p <C I (^ near I, nearly two sepa- 
rate clusters) we find that the curve is well approximated 
by dcx -ln(j3)/(2u- I). 
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SUPPORTING FIG. S6. Variation in the behavior of systems 
(2) and (4) with increasing perturbation off of all-to-all. This 
illustrates delay time as inter-cluster connection probability 
p varies. Equivalent values of the perturbation parameter S 
in order of decreasing p are S = 0, S = 0.14, S = 0.60, and 
S = 0.98. Left panel: Discrete stochastic system (2) (ensemble 
averages of 10 realizations). Right panel: Continuous deter- 
ministic system (4). For all simulations x(0) = 0.1, u = 0.6 
and N = 500. 



7. TIME DELAY 

We define the effective time delay d to be the delay 
between the perturbed solution (not all-to-all) and the 
all-to-all solution (the logistic function, when Pyx — cx'^u 
with a « I), as measured when R has risen halfway to its 
asymptotic value of I (we assume a rising function with 
no loss of generality: the symmetric case of a decaying 
function can be examined under the change of variables 
u i-> I — u, 1-7- 1 — Xq). Thus d is the difference in the 



We find this form by assuming R{t) ^ Rq + rjy{t) + 
0{if') for rj <^1, then eliminating terms of order higher 
than linear in the equation governing R ((4) after simpli- 
fying for two cliques) . We then expand this approximate 
equation for small p, retaining only lowest order terms. 
The resulting system is linear and can be solved exactly 
for the critical time tc at which R rises to (I + xo)/2. 
The delay is simply the difference between that time and 
= — ln(a:o/(l + xq)) / {2u — I), the critical time for the 
all-to-all p = 1 system. 
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Based on numerical work, the general behavior of this 
approximation — d oc — Inp — seems to remain valid even 
for large p, although the additive constant seems to 
change (see Figure S7). That is to be expected, since 
d — > is required for p — > 1 . 
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SUPPORTING FIG. S7. Variation of time delay d with in- 
creasing inter-cluster connection probability p. Points and 
error bars indicate mean and standard deviation with 10 real- 
izations. Lines show estimated logarithmic dependence. Left 
Panel: Discrete stochastic system (2). Right panel: Contin- 
uous deterministic system (4). For all simulations A*' = 500, 
a;(0) = 0.1 and u = 0.6. 



